### This file runs the actual analyses 

rm(list=ls())
setwd("~/dropbox/CET/our products/work in progress/doer til doer/bjps/replication files")
## read in packages 
# install.packages("ri") ## uncomment to install packageslibrary(ri)
library(ri)
# source in data 
source("bhatti_data_recreate.r")

##Set number of iterations for randomization inference
iter <- 100000
set.seed(061015)

##Experiment 1, Copenhagen 

cph_Z <- cph_data$assigned
cph_O <- cph_data$delivered
cph_Y <- cph_data$voted

##turnout for control group 
mean(cph_Y[cph_Z==0])

##Define treatment probabilities
probs <- cph_data$probs

##Define approximate permutation matrix
perms_cph <- genperms(cph_Z, blockvar = probs, maxiter = iter)

##Estimate ITT 

itt_cph   <- estate(cph_Y, cph_Z, prob = probs)

##gen hypothesized PO's
Ys_cph    <- genouts(cph_Y, cph_Z, ate = itt_cph)

##generate randomization distribution of itt's
distout_cph <- gendist(Ys_cph, perms_cph, prob = probs)

##Scale standard error

se_itt_cph <- 
  rep_tabs[[3]] * dispdist(distout_cph, ate = itt_cph)$sd

##Find contact rate
contact_cph <- mean(cph_O[cph_Z==1])

## Find turnout in control group and treatment group
ct_turnout_cph  <- mean(cph_Y[cph_Z == 0])
tr_turnout_cph  <- mean(cph_Y[cph_Z == 1])

# find treatment and control group sizes 
n_cph <- table(cph_Z)

##Summarize randomization distrubution for ITT but dividing by contact rate to obtain CACE-estimate

CACE_cph    <- itt_cph    / contact_cph
CACE_cph_se <- se_itt_cph / contact_cph


# repeat for Randers ------------------------------------------------------
##Experiment 2, Randers 

ran_Z <- ran_data$assigned
ran_O <- ran_data$delivered
ran_Y <- ran_data$voted

##turnout for control group 
mean(ran_Y[ran_Z==0])

##Define treatment probabilities
probs <- ran_data$probs

##Define approximate permutation matrix
perms_ran <- genperms(ran_Z, blockvar = probs, maxiter = iter)

##Estimate ITT 

itt_ran   <- estate(ran_Y, ran_Z, prob = probs)

##gen hypothesized PO's
Ys_ran    <- genouts(ran_Y, ran_Z, ate = itt_ran)

##generate randomization distribution of itt's
distout_ran <- gendist(Ys_ran, perms_ran, prob = probs)

##Scale standard error

se_itt_ran <- 
  rep_tabs[[4]] * dispdist(distout_ran, ate = itt_ran)$sd

##Find contact rate
contact_ran <- mean(ran_O[ran_Z==1])

## Find turnout in control group and treatment group
ct_turnout_ran  <- mean(ran_Y[ran_Z == 0])
tr_turnout_ran  <- mean(ran_Y[ran_Z == 1])

# find treatment and control group sizes 
n_ran <- table(ran_Z)

# find treatment and control group sizes 
table(ran_Z)

##Summarize randomization distrubution for ITT but dividing by contact rate to obtain CACE-estimate

CACE_ran    <- itt_ran    / contact_ran
CACE_ran_se <- se_itt_ran / contact_ran


